Read in data

fert.data <- read_csv(here::here("fert-data.csv"))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   pH.experim = col_double(),
##   Perc.Fertilization = col_double(),
##   Insemination.mins = col_double(),
##   Fert.success.mins = col_double(),
##   Sperm.pre.exp.time = col_double(),
##   egg.pre.exp.time = col_double(),
##   pH.control = col_double(),
##   pH.delta = col_double(),
##   Sperm.per.uL = col_double(),
##   Sperm.per.mL = col_number(),
##   n.females = col_double(),
##   n.males = col_double(),
##   Yeear = col_double()
## )
## See spec(...) for full column specifications.
fert.data.2 <-
  fert.data %>%
  dplyr::select(pH.experim, Perc.Fertilization, 
         Insemination.mins, Fert.success.mins, Sperm.pre.exp.time, 
         egg.pre.exp.time, pH.delta, Sperm.per.mL, sperm.egg, n.females, n.males) %>%
  mutate_if(is.factor, as.numeric) %>%
  mutate_if(is.character, as.numeric)
## Warning: NAs introduced by coercion
fert.data.3 <-
  fert.data %>%
  mutate_at(c("Phylum", "Common name", "Brooders/Spawniers", "Family", "Taxa", "Species") , as.factor)
fert.data.3$pH.group <- cut(fert.data.3$pH.experim,c(6,7.3,7.5,7.8, 8.2))
fert.data.3$Phylum <- factor(fert.data.3$Phylum, levels = c("Echinodermata", "Mollusca","Cnidaria","Crustacean"))
str(fert.data.3)
## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 465 obs. of  26 variables:
##  $ Phylum            : Factor w/ 4 levels "Echinodermata",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Common name       : Factor w/ 31 levels "amphipod","Antarctic bivalve",..: 10 10 10 10 19 19 23 23 21 21 ...
##  $ Brooders/Spawniers: Factor w/ 3 levels "Brooder","Spawner",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Family            : Factor w/ 27 levels "Acartia","Acroporidae",..: 2 2 2 2 18 18 18 18 18 18 ...
##  $ Taxa              : Factor w/ 13 levels "abalone","amphipod",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ Species           : Factor w/ 42 levels "Acartia erythraea",..: 4 4 4 4 16 16 36 36 33 33 ...
##  $ pH.experim        : num  8.15 7.85 8.15 7.85 8.15 7.85 8.15 7.85 8.03 NA ...
##  $ Perc.Fertilization: num  99 90 42 48 87 86 97 98 57 43 ...
##  $ Insemination.mins : num  180 180 180 180 180 180 180 180 420 420 ...
##  $ Fert.success.mins : num  180 180 180 180 180 180 180 180 420 420 ...
##  $ Sperm.pre.exp.time: num  0 0 0 0 0 0 0 0 NA NA ...
##  $ egg.pre.exp.time  : num  0 0 0 0 0 0 0 0 NA NA ...
##  $ pH.control        : num  8.15 8.15 8.15 8.15 8.15 8.15 8.15 8.15 8.03 8.03 ...
##  $ pH.delta          : num  0 -0.3 0 -0.3 0 -0.3 0 -0.3 0 NA ...
##  $ Sperm.per.uL      : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Sperm.per.mL      : num  1e+05 1e+05 1e+05 1e+05 1e+05 1e+05 1e+05 1e+05 3e+06 3e+06 ...
##  $ sperm.egg         : chr  "50000" "50000" "50000" "50000" ...
##  $ n.females         : num  1 1 1 1 1 1 1 1 3 NA ...
##  $ n.males           : num  1 1 1 1 1 1 1 1 3 NA ...
##  $ Other             : chr  "used engauge digitizer" "used engauge digitizer" "used engauge digitizer" "used engauge digitizer" ...
##  $ Data source       : chr  NA NA NA NA ...
##  $ Notes/other treat : chr  "temp = 27" "temp = 27" "temp = 27" "temp = 27" ...
##  $ DOI               : chr  NA NA NA NA ...
##  $ Author            : chr  "M. Schutter, Y. Nozawa, H. Kurihara, 2009" "M. Schutter, Y. Nozawa, H. Kurihara, 2009" "M. Schutter, Y. Nozawa, H. Kurihara, 2009" "M. Schutter, Y. Nozawa, H. Kurihara, 2009" ...
##  $ Yeear             : num  2015 2015 2015 2015 2015 ...
##  $ pH.group          : Factor w/ 4 levels "(6,7.3]","(7.3,7.5]",..: 4 4 4 4 4 4 4 4 4 NA ...
fert.data.4 <- fert.data.3[,c(1,5,7:9,11,12,16:19)] %>%
    mutate_if(is.factor, as.numeric) %>% #convert factors to numeric factors 
    mutate_if(is.character, as.numeric)  #convert character column to numeric 
## Warning: NAs introduced by coercion

Generate distance matrixusing gowers coefficient

Gowers allows for missing data and multiple data types

dist.gower <- vegdist(fert.data.4, "gower") 

Perform PCoA

Its usage is: cmdscale(d, k, eig = FALSE, add = FALSE) where: • d is a dissimilarity object (generated by dist or vegdist) • k is the number of principal components (PC) that should be extracted from the distance matrix (max number = min(col, rows)-1) • eig, logical. If TRUE eigenvalues for each PC are retuned. Default: FALSE. • add, logical. If TRUE a constant is added to each value in the dissimilarity matrix so that the resulting eigenvalues are non-negative. Default: FALSE.

The principal scores are contained in spe.pcoa$points and the eigenvalues are contained in spe.pcoa$eig.

spe.pcoa <- cmdscale(dist.gower, eig=T, add=T, k=2)
head(spe.pcoa$points, n=15)
##             [,1]       [,2]
##  [1,] -0.5743930 -1.0530913
##  [2,] -0.6098390 -0.8683141
##  [3,] -0.7893031 -0.2945215
##  [4,] -0.7845755 -0.2677329
##  [5,] -0.6166947 -0.9347762
##  [6,] -0.6298105 -0.8185384
##  [7,] -0.5788613 -1.0403509
##  [8,] -0.5844064 -0.9437828
##  [9,] -1.0266021 -0.3297779
## [10,] -1.6639959  0.2962041
## [11,] -1.0892730  0.0416065
## [12,] -1.8912539  0.5178418
## [13,] -1.9568934  0.9443114
## [14,] -1.9565642  0.9317737
## [15,] -0.8084025 -0.9088631
head(spe.pcoa$eig, n=15)
##  [1] 327.84640 193.43571 136.24688  87.00823  75.00325  68.77631  67.36747
##  [8]  61.84909  58.84786  53.01254  51.42648  48.04094  44.58181  43.14809
## [15]  40.83116

Calculate the percent of variation explained by principal coordinates:

hist(spe.pcoa$eig/sum(spe.pcoa$eig)*100)

Compare the eigenvalues to expectations according to the broken stick model.

plot(spe.pcoa$eig[1:100]/sum(spe.pcoa$eig)*100,type="b",lwd=2,col="blue",xlab= "Principal Component from PCoA", ylab="% variation explained", main="% variation explained by PCoA (blue) vs. random expectation (red)")
lines(bstick(100)*100,type="b",lwd=2,col="red")

View the ordination plot.

This plot represents each of the sites in 2-D ordination space (x-axis = principal component 1, y-axis = principal component 2). (Should try to use something relevant for text)

Calculate the PC loadings (i.e., variable weights)

Calculate and depict species loadings (i.e., principal weights in the eigenvectors) on each principal coordinate. Use the function envfit() along with the PC scores from our PCoA object. The function envfit() performs a linear correlation analysis based on standardized data (in other words, a simple linear regression) between each of the original descriptors (i.e., species) and the scores from each principal component. A permutation test is used to assess statistical significance, rather than using the F distribution.

print(vec.sp<-envfit(spe.pcoa$points, k=45, fert.data.4, perm=1000, na.rm=T))
## 
## ***VECTORS
## 
##                        Dim1     Dim2     r2   Pr(>r)    
## Phylum             -0.88963 -0.45669 0.8771 0.000999 ***
## Taxa                0.90418  0.42715 0.8888 0.000999 ***
## pH.experim          0.03946 -0.99922 0.2462 0.000999 ***
## Perc.Fertilization  0.22207 -0.97503 0.7852 0.000999 ***
## Insemination.mins  -0.21279  0.97710 0.0960 0.002997 ** 
## Sperm.pre.exp.time -0.83941  0.54349 0.0583 0.041958 *  
## egg.pre.exp.time   -0.78719  0.61671 0.5652 0.000999 ***
## Sperm.per.mL       -0.79476  0.60692 0.5732 0.000999 ***
## sperm.egg          -0.99932  0.03686 0.6199 0.000999 ***
## n.females          -0.66979  0.74255 0.5188 0.000999 ***
## n.males            -0.68888  0.72487 0.2138 0.000999 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1000
## 
## 340 observations deleted due to missingness

Plot the eigenvectors on the ordination plot

p.max is the significance level that the species occurrence data must have with either PC in order to be depicted (these p-values were presented in vec.sp).

fert.data.4$pH.group <- cut(fert.data.4$pH.experim, c(6,7.3,7.5,7.8, 8.2))
pdf(file = "fert.PCoA.pdf", width = 9, height = 8)
pl <- ordiplot(spe.pcoa, type = "none",  xlim = c(-1,1.5))
## species scores not available
points(pl, "sites", cex=0.8, pch=c(21,22,23,24)[fert.data.4$pH.group], 
       bg=c("red","blue","green","purple")[fert.data.4$Phylum])
plot(vec.sp, p.max=.01, col="black") # note, p.max = set p-value threshold for plotting vectors. 
legend(x="topright", legend = levels(fert.data.3$Phylum), col=c("red","blue","green", "purple"), pch=c(16,16,16,16))
legend(x="right", legend = levels(fert.data.3$pH.group),pch=c(21,22,23,24))
dev.off()
## quartz_off_screen 
##                 2

Re-run PCoA, leave pH out of matrix, but then color code points that way.

dist.gower2 <- vegdist(fert.data.4[c(2,4:11)], "gower", na.rm = F) 
spe.pcoa2 <- cmdscale(dist.gower2, k=2, eig=T, add=T)
print(vec.sp2<-envfit(spe.pcoa2$points, k=45, fert.data.4[c(2,4:11)], perm=1000, na.rm=T))
## 
## ***VECTORS
## 
##                        Dim1     Dim2     r2   Pr(>r)    
## Taxa                0.63598 -0.77170 0.8539 0.000999 ***
## Perc.Fertilization  0.61708  0.78690 0.9546 0.000999 ***
## Insemination.mins  -0.85462 -0.51925 0.0472 0.054945 .  
## Sperm.pre.exp.time -0.98964  0.14356 0.0488 0.065934 .  
## egg.pre.exp.time   -0.99754  0.07008 0.5951 0.000999 ***
## Sperm.per.mL       -0.99696  0.07787 0.6043 0.000999 ***
## sperm.egg          -0.92696  0.37515 0.5829 0.000999 ***
## n.females          -0.99770 -0.06773 0.5179 0.000999 ***
## n.males            -0.99765 -0.06853 0.1923 0.000999 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 1000
## 
## 340 observations deleted due to missingness
pdf(file = "fert.PCoA-nopH.pdf", width = 9, height = 8)
pl2 <- ordiplot(spe.pcoa2, type = "none",  xlim = c(-1,1.5))
## species scores not available
points(pl2, "sites", cex=0.8, pch=c(21,22,23,24)[fert.data.4$pH.group], 
       bg=c("red","blue","green","purple")[fert.data.4$Phylum])
plot(vec.sp2, p.max=.01, col="black") # note, p.max = set p-value threshold for plotting vectors. 
legend(x="topright", legend = levels(fert.data.3$Phylum), col=c("red","blue","green", "purple"), pch=c(16,16,16,16))
legend(x="right", legend = levels(fert.data.3$pH.group),pch=c(21,22,23,24))
dev.off()
## quartz_off_screen 
##                 2

Run correlation test with factors with significant loadings on dimensions 1 and 2

All of the below except Sperm.pre.exp.time, sperm.egg, and n.males.

VECTORS Dim1 Dim2 r2 Pr(>r)
pH.experim -0.78537 0.61903 0.9157 0.000999
Perc.Fertilization -0.76812 -0.64030 0.8868 0.000999 Insemination.mins 0.42319 -0.90604 0.3439 0.000999 Fert.success.mins 0.38683 -0.92215 0.3583 0.000999 Sperm.pre.exp.time 0.95118 0.30863 0.0034 0.880120
egg.pre.exp.time 0.82583 0.56392 0.2503 0.001998
pH.delta -0.73572 0.67729 0.7674 0.000999
Sperm.per.mL 0.81881 0.57406 0.2492 0.001998 sperm.egg 0.20294 0.97919 0.0953 0.038961 *
n.females 0.90918 0.41640 0.2463 0.001998 ** n.males 0.86065 0.50920 0.0570 0.127872

Check out correlations among variables

  • Insemination minutes ~ Fertilization success mins - only use insemination minutes
  • pH delta ~ pH experimental - only use pH experimental
    other correlations, but not relevant (e.g. sperm/mL and egg pre-exposure time)
pdf(file = "fert-correlation-panel.pdf", width = 12, height = 8.5)
pairs(na.omit(fert.data.4), lower.panel=panel.smooth, upper.panel=panel.cor) 
dev.off()
## quartz_off_screen 
##                 2

Plot factors that were deemed sign. in PCoA

Test experimental pH

pH.experim -0.78537 0.61903 0.9157 0.000999 ***

# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=pH.experim, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
  ggtitle("Fertilization Rate ~ pH exposure by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 6 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*pH.experim, fert.data.3)) #significant alone and interaction  
##                            Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Phylum)              3  15799    5266   7.590 5.82e-05 ***
## pH.experim                  1  76312   76312 109.982  < 2e-16 ***
## factor(Phylum):pH.experim   3  10300    3433   4.948  0.00217 ** 
## Residuals                 451 312932     694                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 6 observations deleted due to missingness

Test insemination minutes

Insemination.mins 0.42319 -0.90604 0.3439 0.000999 ***

# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=Insemination.mins, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
  ggtitle("Fertilization Rate ~ Insemination minutes by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 65 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*Insemination.mins, fert.data.3)) # not sign. alone, but sign for some phyla 
##                                   Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Phylum)                     2  22015   11007  14.131 1.18e-06 ***
## Insemination.mins                  1      6       6   0.008   0.9288    
## factor(Phylum):Insemination.mins   2   6310    3155   4.050   0.0182 *  
## Residuals                        394 306910     779                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 65 observations deleted due to missingness

Test egg pre-exposure time

egg.pre.exp.time 0.82583 0.56392 0.2503 0.001998 **

# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=log(egg.pre.exp.time), y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
  ggtitle("Fertilization Rate ~ Egg pre-exposure minutes by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 266 rows containing non-finite values (stat_smooth).
fert.data.3$egg.pre.exp.time.log <- log(fert.data.3$egg.pre.exp.time+1)
summary(fert.data.3$egg.pre.exp.time.log)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   2.773   1.991   2.773   7.273     156
summary(aov(Perc.Fertilization ~ factor(Phylum)*egg.pre.exp.time.log, fert.data.3)) # not sign. alone, but sign for some phyla
##                                      Df Sum Sq Mean Sq F value   Pr(>F)
## factor(Phylum)                        2   2712    1356   1.911 0.149717
## egg.pre.exp.time.log                  1    528     528   0.744 0.389184
## factor(Phylum):egg.pre.exp.time.log   1   8413    8413  11.855 0.000656
## Residuals                           304 215740     710                 
##                                        
## factor(Phylum)                         
## egg.pre.exp.time.log                   
## factor(Phylum):egg.pre.exp.time.log ***
## Residuals                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 156 observations deleted due to missingness

Test sperm concentration

Sperm.per.mL 0.81881 0.57406 0.2492 0.001998 **

ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=log(Sperm.per.mL+1), y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
  ggtitle("Fertilization Rate ~ sperm concentration (log-trans) by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 204 rows containing non-finite values (stat_smooth).
fert.data.3$Sperm.per.mL.log <- log(fert.data.3$Sperm.per.mL+1)
summary(aov(Perc.Fertilization ~ factor(Phylum)*Sperm.per.mL.log, fert.data.3))  # not sign. alone, but sign for some phyla
##                                  Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Phylum)                    2  29032   14516  27.612 1.40e-11 ***
## Sperm.per.mL.log                  1    537     537   1.022    0.313    
## factor(Phylum):Sperm.per.mL.log   2  23492   11746  22.343 1.15e-09 ***
## Residuals                       255 134059     526                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 204 observations deleted due to missingness

Same plot as above, but shapes = pH group

pH groups: – 6-7.3 – 7.3-7.5 – 7.5-7.8 – 7.8-8.2

ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=log(Sperm.per.mL+1), y=Perc.Fertilization, group=Phylum:pH.group, col=Phylum:pH.group, text=`Common name`, shape=pH.group)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum:pH.group)) +
  ggtitle("Fertilization Rate ~ sperm concentration (log-trans) by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 204 rows containing non-finite values (stat_smooth).
## Warning: Factor `Phylum:pH.group` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Warning: Factor `Phylum:pH.group` contains implicit NA, consider using
## `forcats::fct_explicit_na`

Test sperm: ratio

sperm.egg 0.20294 0.97919 0.0953 0.038961 *

# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=sperm.egg, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
  ggtitle("Fertilization Rate ~ sperm:egg ratio by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 2 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*sperm.egg, fert.data.3))  # sign. main and interaction effects 
##                           Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Phylum)             2  16503    8252   29.06 1.16e-11 ***
## sperm.egg                 35 135929    3884   13.68  < 2e-16 ***
## factor(Phylum):sperm.egg   1   9519    9519   33.52 3.09e-08 ***
## Residuals                180  51119     284                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 246 observations deleted due to missingness

Test number of females used in assays

n.females 0.90918 0.41640 0.2463 0.001998 **

# plot % fert ~ pH.experim by Phylum
ggplotly(fert.data.3 %>%
ggplot(mapping=aes(x=n.females, y=Perc.Fertilization, group=Phylum, col=Phylum, text=`Common name`)) + 
  geom_point(size=1.5, width=0.02) +
  #facet_wrap(~Taxa) +
  geom_smooth(method="lm", se=TRUE, aes(fill=Phylum)) +
  ggtitle("Fertilization Rate ~ No. females used for assay, by phylum"))
## Warning: Ignoring unknown parameters: width
## Warning: Removed 97 rows containing non-finite values (stat_smooth).
summary(aov(Perc.Fertilization ~ factor(Phylum)*n.females, fert.data.3))  # sign. main effect, not interaction  
##                           Df Sum Sq Mean Sq F value   Pr(>F)    
## factor(Phylum)             2  14574    7287   9.605 8.62e-05 ***
## n.females                  1  12866   12866  16.957 4.74e-05 ***
## factor(Phylum):n.females   2   2094    1047   1.380    0.253    
## Residuals                362 274654     759                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 97 observations deleted due to missingness

Full model, all factors explored above.

summary(aov(Perc.Fertilization ~ factor(Phylum)*(pH.experim + Insemination.mins + egg.pre.exp.time.log + Sperm.per.mL.log + n.females), fert.data.3))  #  
##                                      Df Sum Sq Mean Sq F value   Pr(>F)
## factor(Phylum)                        2   3959    1979   7.094 0.001069
## pH.experim                            1   3416    3416  12.242 0.000582
## Insemination.mins                     1    153     153   0.548 0.459852
## egg.pre.exp.time.log                  1   1474    1474   5.284 0.022615
## Sperm.per.mL.log                      1   4207    4207  15.079 0.000142
## n.females                             1   7482    7482  26.815 5.68e-07
## factor(Phylum):pH.experim             2    995     497   1.782 0.171047
## factor(Phylum):Insemination.mins      1  10689   10689  38.308 3.64e-09
## factor(Phylum):egg.pre.exp.time.log   1   2883    2883  10.333 0.001536
## factor(Phylum):Sperm.per.mL.log       1   3397    3397  12.173 0.000603
## factor(Phylum):n.females              1    794     794   2.847 0.093164
## Residuals                           190  53015     279                 
##                                        
## factor(Phylum)                      ** 
## pH.experim                          ***
## Insemination.mins                      
## egg.pre.exp.time.log                *  
## Sperm.per.mL.log                    ***
## n.females                           ***
## factor(Phylum):pH.experim              
## factor(Phylum):Insemination.mins    ***
## factor(Phylum):egg.pre.exp.time.log ** 
## factor(Phylum):Sperm.per.mL.log     ***
## factor(Phylum):n.females            .  
## Residuals                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 261 observations deleted due to missingness